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^ ' Abstract 
'■ 

Analytical solutions of the two dimensional triangular and square lattice Boltzmann 
. BGK models have been obtained for the plain Poiseuille flow and the plain Couette flow. 

The analytical solutions are written in terms of the characteristic velocity of the flow, 
the single relaxation time r and the lattice spacing. The analytic solutions are the exact 
. representation of these two flows without any approximation. 

o , 

5 ■ 1 Introduction 
O 

ly^ ■ Since the appearance of lattice gas automata (LGA) and its later derivative, the lattice 

. Boltzmann equation method (LBE) as alternative computational methods to study trans- 

^ I port phenomena, some analytical solutions have been obtained for these methods for nonuni- 
form flows in 2D or 3D models |l|, |2|, ^, ^. They are based on linearized Boltzmann models, 



and they have employed approximations. In IQ, ^, the solution around a global equi- 
librium with constant density and isotropic velocity (zero velocity) was considered. In Q, 
the first order and second order of deviation of distribution function from equilibrium were 
O ■ assumed to take a certain form in terms of flow quantities, and the coefficients in this form 

were obtained using a Chapman-Enskog procedure. These analytic results provided insight 
for applications of the methods. For example, the analytical results allow one to calculate 



5_j ■ viscosity from given collision rules and to estimate and to improve bounce-back boundary 



conditions. They are valuable in enhancing our understanding of the method. Neverthe- 
less, analytical solutions for real flows with boundaries like the Poiseuille flow, which is 
represented exactly by a second-order finite-difference scheme on a uniform mesh, have not 
been obtained previously for LGA or LBE. One reason may be that the boundary condi- 
tions used in LGA and LBE are not exact. For example, bounce-back or a combination of 
bounce-back and specular reflection for modeling the non-slip boundary condition are 
only approximate. The effective non-slip boundary is inside the bounce-back row |2|, |^. 
Recently, Noble et al. have proposed a boundary condition for the lattice Boltzmann 
BGK model (LBGK) on a triangular lattice. When this boundary condition was applied to 
plain Poiseuille flow, the steady state solution of the distribution function gave a parabolic 
velocity profile up to machine accuracy. The result suggests the existence of an analytical 
solution to LBGK, which is an exact representation of the Poiseuille flow. This analytical 
solution is derived in this paper together with an analytical solution for the plain Couette 
flow. 
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2 Analytical Solutions of the triangular lattice LBGK model 



First let us consider the lattice Boltzmann model on a triangular lattice. For a channel flow, 
a triangular lattice is constructed as shown in Fig. 1. There are two types of particles on 
each node of an FHP model: rest particles (type 0) with eg = and moving particles (type 
1 to 6) with unit velocity = (cos((z — l)7r/6), sin((z — l)7r/6), i = 1, ■ ■ ■ ,6 along 6 direc- 
tions. Consider the particle distribution functions fi{x,t), which represent the probability 
of finding a particle at node x = (x,y) and time t with velocity Cj. The lattice Boltzmann 
BGK model proposed in |7| is the equation for the evolution of ff. 

/.(x + 5e„t + 5)-/,(x,t) = -i[/i(x,i)-/f (x,t)], i = 0,---,6 (1) 

where fji'^\x, t) is the equilibrium distribution of the particle of type i at x, t, the right hand 
side represents the collision term and r is the single relaxation time which controls the rate 
of approach to equilibrium. The density per node, p, and the macroscopic flow velocity, u, 
are defined in terms of the particle distribution function by 

6 6 

i=0 i=l 

The equilibrium distribution functions depend only on local density and velocity. A suitable 
equilibrium distribution for the FHP model can be chosen in the following form [^: 

Jo = do - pu =ap- pu , 

1 1 

=d+-p[(ei-u) + 2(ei-u)2--u-u], i = l,---,6, (3) 

where a is an adjustable parameter, d = and Yl,i fi^'^ = P^ Z]j fi^^^i = P^- Note that 

Eq. (H) is written in physical units with the value of the lattice link being 6, Using unit 
speed for particles (with some physical time unit), a time step has the value of 6 as well. A 
Chapman-Enskog procedure can be applied to Eq. (|l]) to derive the macroscopic equations 
of the model. They are given by: the continuity equation (with an error term 0{6'^) being 
omitted): 

^ + V-{pu)=0, (4) 
and the momentum equation (with terms of 0(5^) and 0{6u^) being omitted): 

dtipUa) + di3{pUaUp) = -da{clp) + dfi{2v pSap) , (5) 

where the Einstein summation convention is used, S^p = \[daUp + dpUa) is the strain-rate 

tensor, the pressure is given by p = c^/>, where Cs is the speed of sound with = — ^ — , 
2t - 1 

and V = 5, with v being the the kinematic viscosity. The form of the error terms 

8 

and the derivation of these equations can be found in |Q, |9| . The macroscopic equations of 
LBGK represent the incompressible Navier-Stokes equations in the limit as (5 — > 0, p — > po 
(a constant ) and the Mach number approaches zero. 
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The plain Poiseuille flow in a channel with width 2L and velocity u = (u, v) is given by: 

„ = ,,„(1-|!), „ = 0, | = -0, 1 = 0. (6) 

where G is a constant related to the characteristic velocity uq by 

G = 2pI^no/L^ (7) 

and the flow density p is a constant. This is an exact solution of the incompressible Navier- 
Stokes equations: 

V • u = 0, 

dtipUa) + disipUaUp) = -d^P + H-dpisUa, (8) 

where ^ = pu. Without loss of generality, we assume that L = 1 (a simple scaling of 
y' = y/L makes y' £ [—1,1]). To approximate Poiseuille flow using the lattice Boltzmann 
model, it is convenient to replace the constant gradient by a body force g so that = —Vp. 
The momentum equation of N-S equations with a body force is: 

dtipUa) + dp{pUaUi3) = -daP + PQa + pdpfjUa- (9) 

The Poiseuille flow can be generated with a body force g with pg^ = G,gy = 0, where the 
pressure is held constant. An LBGK which incorporates the body force is a modification of 
Eq. (|^) given by: 

fi{^ + 6ei,t + 5) - fi{^,t) = -^[fi{^,t)-f\^,t)]+5hi, (10) 

r 

where the hi are chosen as: 

ho = 0, hi = +G/A, i = 1,2,6, hi = -G/4, i = 3,4,5, (11) 

so that 



^hi = 0, hiGi = pg, hiCiaCifs = 0. 



Now suppose there exists a solution fi{x,t) of Eq. ( pUj ) and it exactly represents the 
Poiseuille flow. We then expect the following properties: 

1. /j(x, t) is steady (independent of t). 

2. fi{x,t) is independent of x, hence is only a function of y, denoted by fi{y). 

3. /2(y) = fei-y), and fsiy) = f^i-y) from the symmetry of the flow. 

4- T,ifi{y) = P (constant). 

5- Ei fi{y)eix = pu{y) , where u{y) = uq{1 - y^) (remember L = 1). 

6- T,ifi{y)eiy = . 

According to Eq. (^, the equilibrium distributions are given by 

fo^=do-pv?, {u = u{y)) 
ff' = d-^^u, ff^=d+P^u. (12) 
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Using properties 1,2 and Eq. (|10D for i = gives: 

/o(y) = /o(y)-i(/o(y)-/r(2/)), 

hence 

fo{y) = f!>'\y) = do-pu^. (13) 

Eq. (10) for i = 1 gives: 

fiiy) = fi{y) - ^(/i(y) - A'\y)) + SG/i, 

hence 



Similarly 



My) = fi'\y)+T6G/A. (14) 

My) = fPiy)-r5G/4. (15) 

It is seen that /o, /i, /2 are functions of y^, through dependence of u and u^. To find the 
remaining fi{y), we note that ff^\i = 2,3,5,6 have no term and thus are functions of 
only, so the following form is suggested: 

fi{y) = ai + biy + Ciy'^, i = 2,3, 5,6, (16) 

where the twelve unknown coefficients Oj, bi,Ci depend on flow quantities, r, dy, but not on 
y. Using property 3, we obtain 

02 + b2y + C22/^ = ae - bey + cqv^, as + bsy + cs?/^ = 05 - 65?/ + C5y^, 
which should be true for any y, thus: 

06 = 02, &6 = -b2, C6 = C2 

05 = 03, ^5 = -^3, C5 = C3. (17) 

Similarly, using property 6, we find: 

2b2y + 2b^y = 0, hence 63 = -62. (18) 
Property 4 with information obtained gives: 

2(02 + 03) + 2(c2 + C3)y2 + Uy) + h{y) + U{y) = P. (19) 
On using the expressions for /o, /i, /4 given in Eq. (|13| - |15|), we find 

02 + 03 = 2 (/^ ~ '^0 ~ ~ (^'^ using the expression of d), C3 + C2 = 0; (20) 
which gives 

C3 = -C2, a3 = 2d-a2. (21) 
Similarly, property 5 with information obtained yields: 

a2 = l:puo + d-TSG/4; C2 = -\puQ. (22) 
6 6 
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At this point, only 62 remains unknown. Using Eq. ( |10D for z = 2, we have 

/2(y + dy) = h{y) - i(/2(y) - f^2\y)) + (23) 

where dy is the vertical spacing between two lattice rows and dy = (\/3/2)5. O n using the 
expression for 72''^ we can obtain: 

1 11 

C2[y'^ + 2ydy + {dyf\ + b2y + h2dy + a2 = (1 - -)(c2y^ + 62^ + 02) + -/9Uo(l -y^)] +5G/4. 

r r b 

(24) 

The balance of terms linear in y yields: 

62 = -2TC2dy = ^Tpuody, (25) 

and fortunately the equations for the coefficients of and y^ are both satisfied. It is easy 
to check that the evolution equations for fs, f^, fe are all satisfied with the choice of ai,bi, Ci 
obtained so far. 

Putting these results all together, we find that the quantities 

fo = do- pu^, 

f, = d+^u+^u^ + T6G/4, 

h = d-^u+^u^-r5G/i, 

11 1 

/2 = --^puoy^ + -Tpuoydy + -puo + d- t6G/4, 

DO O 
11 1 

h = +-^puoy^ - -Tpuoydy - -puo + d + t6G/4, 

DO O 
11 1 

h = +-^puoy'^ + -Tpuoydy - -puo + d + t6G/4, 

DO O 

fe = --^PUQV^ - i^rpuoydy + -puo + d- t5G/A, (26) 

DO D 

all satisfy properties 1-6 and that they together with the equilibrium distribution given in 
Eq. (12) satisfy the LBGK equation, Eq. (^). Hence it is an exact representation of the 
Poiseuille flow. 

Next, let us see to what boundary condition the solution in Eq. (^) corresponds. Taking 
the bottom boundary with y = —1, u = 0, we have 

/o = do, h = d + T5G/A, U = d-T5G/A 
f2 = --Tpuody + d- t5G/A, /s = +-Tpuody + d + t5G/A 

h = -^Tpuody + d + t6G/4, h = +^Tpuody + d- t6G/4. (27) 

It is seen that on the bottom, after the collision and forcing, /2 = /s — 2t5G/4:, = 
/g + 2t5G/4. Hence, if a bounce-back boundary condition /2 = /sj/s = /e, on /2,/3 is 
applied at the bottom to replace the collision and forcing step, the error is of order 5. 
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This shows that the bounce-back boundary condition is first-order accurate. This has been 
confirmed in computations |l^, |ll| . 

To obtain the steady-state analytical solution in the LBGK simulations, the boundary 
condition should be suitably chosen for the simulation. The boundary condition proposed 
by Nobel et al. ||5| is a suitable choice. If we are looking at a node B on the bottom, after 
streaming, /2 and /s are empty at the node B since no particle is coming from outside. 
Then Eq. (|2|) with u = v = are used to determine p,f2,f3- Then the normal collision 



with force as given in Eq. (10) is applied to fi on the boundaries. Suppose that initially, 
we use uniform density po and zero velocity through the flow field; then we compute fj^'^^ 

and set fi = f-^^ through the field. Since there are no pressure gradients, it is natural 
that the density at each node remains constant pQ (confirmed by simulations). Therefore 
Eq. (|2|) can be used to find the unique /2 , /s with the correct density and velocity, hence it is 
consistent with the evolution of /2, /s in the analytical solution. Simulation results indicate 
that the numerical solution with this boundary condition approaches the analytical solution 
as t oo. 

Next, let us consider a plain Couette flow, where the flow between two parallel plates 
(corresponding to y = and y = 1) are driven by the constantly moving top plate with 
velocity uq. In this case, the solution is given by: 

u = uoy, 0<y <1, v = 0, Vp = 0, (28) 

with p a constant and with no body force. So the LBGK model Eq. (|l|) is used. Using a 
similar procedure, we find the analytical solution of Eq. (||) representing the Couette flow: 

f J 2 f J , P , P 2 f ^ P , P 2 

Jo = do- pu , ji = d + -u + -u , f4 = d--u + -u, 
/2 = +lpuoy + d- ^Tpuody, h = - \puQy + d+ ^rpuody, 

DO DO 

fb = --^puoy + d+ -Tpuody, fe = +-puoy + d - -rpuody. (29) 

DO 

We note that these analytical solutions are valid for any uq, t and dy. 



3 Analytical Solutions of the square lattice LBGK model 

The square lattice Boltzmann BGK model is proven more robust than the triangular model 
in numerical simulations ||9|, |l^ , it is important and interesting to find analytical solutions 
for it. The square lattice Boltzmann BGK model uses 3 types of particles. Particles of 
type 1 move along the x axis or the y axis with speed ej = (cos(7r(i — l)/2), sin(7r(i — 
l)/2),i = 1,2,3,4, and particles of type 2 move along the diagonal directions with speed 
Bi = \/2(cos(7r(z - 4 - i)/2),sin(7r(i - 4 - ^)/2),i = 5,6,7,8. Rest particles of type 
with eg = (speed zero) are also allowed at each node. Each node is connected to its 8 
nearest neighbors by 8 links of length S (in physical units) or ^/26 as shown in Fig. 2. The 
single-particle distribution function /j(x, t) again satisfies the LBGK model Eq. (|) (with 
i = 0, • • • ,8). The density p and the macroscopic velocoty u are still defined in Eq. (|2|). 
For the square lattice, the equilibrium distribution can be chosen in the following form for 
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particles of each type (the model d2q9 0): 

/f) = ip[l + 3(e,.u) + ^(ei.u)2-^u.u], i = l,2,3,4 

1 9 3 

= ^P[l + 3(e, • u) + -{ei . u)2 - -u • u], i = 5, 6, 7, 8, (30) 

= pu. The macroscopic equations of the model are the 

o' i cr i 

2r — 1 

same as given in Eqs. (|^,^) with = 1/3, and v = — ^ — 5. To incorporate a body force 

in the model to model Poinseuille flow, Eq. (^) ( i = 0, • • • , 8) is used, with hi chosen in 
the following way |13|: 

/io = 0, ^1 = ^^, /i2 = 0, /i3 = -iG, /i4 = 0, 

h^ = H = ^G, he = h = -^G, (31) 



To derive an analytical solution of Eq. (|1C|) for the square lattce, we note that the six 
properties in Section 2 still applies except that property 3 is replaced by: 

3- /2(y) = h{-y), h{y) = h{-y) and /6(y) = h{-y) from the symmetry of the flow. 

Using a similar procedure as in Section 2, we can find that 

fo{y) = fl>'\y) = Ip{1 - lu'), in = no(l - y')), (32) 
h{y) = ^pil + 3u + 3u^) + ^Tiypuo5. (33) 
fsiy) = ^P(l - 3n + 3u^) - ^Ti^puoS. (34) 

fi{y) = ai + biV + CiV^ + div'^ + Civ'^, i = 2, 4, 5, 6, 7, 8, (35) 



and 

with 



04 = a2, 64 = -62, C4 = C2, (i4 = -^2, 64 = -62, 

as = as, ^8 = -^5, C8 = C5, ds = -ds, eg = -65, 

a^ = gq, bj = -be, cy = cq, dy = -de, ey = -65, (36) 



and 



02 = -4t^/3Uo5'' + GT^puld"^ - ^-t'^puIS'^ + '^t'^ pul5^ + ^rpulS^ - ^rpulS^ ^^^O' 
62 = 4r3pug53 _ 4^2 ^^2^3 ^ ^^p^2^3 _ IrpulS, 

C2 = -2t'^pu15'^ + Tpuld"^ + ^png, d2 = '^rpu^d, 62 = -^png, (37) 
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1 2r2 1 .2 1 1 2 1 1 

-tpuq6 +—Tpuo6 + —p + —puQ + —puo + -Tupuod, 

1 11 

b5 = -2T^pul5^ + 2t'^pu16^ - -rpu^d^ + -rpu^d + -rpuoS, 

C5 = r'^pulS^ - ^rpulS"^ - ^pul - j^puo, ^5 = -^rpul5, 65 = j^pul, (38) 

and 

ae = 2r'pul6' - Sr^pulS' + -r^pulS' - pul5' + \r^pu,6^ - ^rpuld' + 

fa O D IZ 

1 2r2 1 .2 1 1 2 1 1 

-tpuq6 - —rpuod + + J^P'^o - J^P^o ~ g^^/^^o*), 

1 11 

be = -2T^pul5^ + 2T^pul6^ - -rpuld^ + -rpuld - -rpuoS, 

C6 = T^pul6'^ - ^Tpul6'^ - ^pul + Y^puo, da = -^Tpul6, = ^pul- (39) 

Eqs. (p2|,pS J34) together with Eqs. (p5,36,37,3^j39|) completely specify the analytical solu- 
tion, which is a solution of Eq. (|^) and it exactly represents the Poiseuille flow. 

Next, let us see to what boundary condition this analytical solution corresponds. Taking 
the bottom boundary with y = — 1, u = 0, we find the relation of fai after the collision and 
forcing: 

4 2 
h- h = ^riypuod, /2 - /4 = -26^{4T^pul - ir^pul + -Tpu^), 

2 1 2 

h- h = --t'^PUqS'^ + -rpulS'^ + 4.T^pul5^ - Ar^pulS^ + -rpu^d^, 

k-h = +\t^Pu^5^ - \tpu15^ + AT^pulS^ - Ar^pu^S^ + ^rpu^d^ (40) 

if a bounce-back boundary condition in which fi exchanges with /3, f2 = fi, fb = fr, /e = 
fs, is applied at the bottom to replace the collision and forcing step, the error introduced to 
/i and /s is of order 5. This shows that the bounce-back boundary condition is first-order 
accurate. This has been confirmed in computations ||9|, |ll|. To obtain the steady-state 
analytical solution derived in this paper in LBGK simulations, the boundary condition 
should be suitably chosen. No numerical simulation on a square lattice Boltzmann BGK 
model has obtained an exact solution for the Poiseuille flow so far. Of course, if we use the 
analytical solution derived as the initial condition and provide the right form of distribution 
functions on the boundaries, we will be able to obtain the analytical solution in a simulation 
(confirmed in simulations). Specification of the analytical solution on boundary does not 
provide a boundary condition of general purpose. Nevertheless, the analytical solution will 
give some guidance in developing better boundary conditions of general purpose for the 
model. 

The analytical solution for plain Couette fiow are given by 
/o = ^P(l - ^^^ol/^)' /i = ^P(l + 3noy + 3uoy^), /s = ^/)(1 - 3uoy + 3ngy2), 
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/2 = —r'^pulS'^ + ^TpulS'^ ^^^^ ^TpulSy - ^puly^, 

h = —t'^pu^S'^ + ^rpuld"^ ^'^P^o^y ~ ^/^^oy^ 

/s = ^r'^puld'^ - ^rpuld'^ - ^rpuod + + {-^'^pu'^S + ^puo)y + ^pu'^y'^, 
/e = ^T^pu^o^^ - ^Tpul6^ + ^rpuoS + ^P + (-^^P^l^ - ^puo)y + ^puly^, 
fi = \t'^PuIS'^ - ^rpulS"^ - ^rpuoS + ^P + (+^^P«o<^ - ^P^^o)?/ + j^puW, 
fs = ^r'^pulS'^ - ^rpulS^ + ^rpuoS +^P+ i+^Tp^lS + ^puQ)y + ^pu^y^. (41) 

For the Couette flow, the top boundary is a moving one, the analytical solution given here 
will give a guidance in developing a suitable boundary condition for moving boundaries. 

We note that these analytical solutions are valid for any uo, r and 6. They will enhance 
our understanding of the method and will give guidance in applications. 
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Figure captions 

Figure 1. The geometry of the plain channel flow. 
Figure 2. Schematic of a square lattice. 
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